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Abstract 

We investigate the DFT+U approach as a viable solution to describe the low-lying states of 
ligated and unligated iron heme complexes. Besides their central role in organometallic chemistry, 
these compounds represent a paradigmatic case where LDA, GGA, and common hybrid functionals 
fail to reproduce the experimental magnetic splittings. In particular, the imidazole pentacoordi- 
nated heme is incorrectly described as a triplet by all usual DFT flavors. In this study we show 
that a U parameter close to 4 eV leads to spin transitions and molecular geometries in quantitative 
agreement with experiments, and that DFT+U represents an appealing tool in the description 
of iron porphyrin complexes, at a much reduced cost compared to correlated quantum-chemistry 
methods. The possibility of obtaining the U parameter from first-principles is explored through a 
self-consistent linear-response formulation. We find that this approach, which proved to be suc- 
cessful in other iron systems, produces in this case some overestimation with respect to the optimal 
values of U. 
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I. INTRODUCTION 



Enzymatic sites containing transition metals are among the most relevant biophysical 
systems currently studied using first-principles quantum mechanical approaches. The ap- 
plication of such tools, however, is often severely limited as a consequence of the inability 
of conventional electronic structure methods — such as Hartree-Fock or density-functional 
theory — to provide a qualitatively correct description of the spin-state energetics of the 
metal center. Iron porphyrins, which constitute the prosthetic group of the ubiquitous heme 
proteins, are a paradigmatic example where the aforementioned approaches can not be relied 
upon to predict the ground state multiplicity of the system. 

The spin state of iron porphyrins, as much as the spin state of any transition metal 
complex, is determined by the coordination symmetry and the nature of the ligands. The 
three lowest accessible spin states (a singlet, a triplet and a quintuplet if the number of 
electrons is even, or a doublet, a quartet and a sextet if it is odd) are conventionally referred 
to as low, intermediate, and high-spin. In unligated porphyrins the metal is coordinated 
to four in-plane nitrogen atoms, and experimental studies on model compounds, namely on 
Fe(II) tetraphenylporphine (FeTPP), indicate for this coordination mode a triplet ground 
state. Additional axial ligands produce alternative multiplicities: imidazole gives rise to 
high spin hemes, while strong ligand-fields as that of diatomic molecules like CO, NO or 
O2, favor low spin configurations.^ Six coordinated hemes, with two axial ligands, usually 
exhibit a low spin state unless the ligand-field is extremely weak. Fig. [T] depicts a schematic 
view of the d-states energy levels in three distinctive coordination environments. 

Even though first-principles approaches, specifically Hartree-Fock (HF) and density- 
functional theory (DFT), greatly contributed to the interpretation and understanding of 
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FeTPP FeTPP(lnn) FeTPP(lm)(CO) 



FIG. 1: Schematic representation of the d-orbitals energy levels for FeTPP. From left to right: free 
(four coordinated), ligated to imidazole (five coordinated), and ligated to imidazole plus CO (six 
coordinated). 

the functional aspects of the active site of heme proteins at the molecular level, attempts 
to predict the ground state multiplicity of these systems soon made apparent that an ac- 
curate description of the electronic structure might require more sophisticated techniques. 
This fact can be tracked down to the spin-transition energies provided by HF and DFT for 
isolated iron atoms and ions or different iron compounds, where it has been systematically 
observed that HF favors high-spin electronic configurations while DFT exhibits a preference 
for low-spin states.-"— Such biases are similarly manifested in heme complexes: Table I 
summarizes this trend in five and six coordinated iron porphines (FeP). 

For the last decade DFT has been the first method of choice to perform electronic struc- 
ture calculations of biological models, and in particular of heme systems. In this context, 
one of the most crucial failures of common exchange-correlation functionals has been de- 
tected in the deoxygenated active site of hemoglobin and myoglobin (Table I). The earliest 
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study reporting this flaw is due to Rovira et al.,-^ who obtained for the five coordinated 
model Fe^^P(Im) (axial ligand: imidazole) a triplet state 6.5 kcal/mol below the quintuplet, 
which is the experimental ground state of the system. After this work, a few others followed 
which also observed this inversion using B3LYP or different pure GGA functionals.-*^ Liao 
and Scheiner claimed to have found a quintuplet ground state for this five coordinated 
compound employing a DFT-GGA functional. In their calculations, however, electronic 
symmetry constraints were imposed. To the best of our knowledge, DFT functionals yield 
for Fe^^P(lm) a triplet ground state in the absence of symmetry constraints. In an effort 
to quantify the errors in the DFT estimates of spin transition energies, Ghosh and Taylor 
resorted to highly-correlated techniques as CASPT2 and CCSD(T) to explore the iron (III) 
porphyrin chloride.^ This is another example of a high spin five coordinated heme complex 
for which DFT predicts a quartet favored over the sextet, in this case by around 7 kcal/mol. 
B3LYP, on the other hand, finds about the same energy for both spin configurations. The 
more accurate approaches CASPT2 and CCSD(T) agree in yielding a sextet ground state, 
16 kcal/mol below the quartet.® A latter work by these authors shows the same low-spin 
bias in B3LYP for the iron (IV) porphyrin difiuoride.- It is worth noting here that even 
CASPT2 — employed with the moderate active spaces currently affordable — has been found 
fallible in the estimation of these elusive spin states. Inaccuracies have been reported in 
the prediction of the electronic ground states of the isolated iron porphyrin^ and of the 
oxyheme.— 

It is possible to find a rationale for the biases in DFT and HF, considering the bal- 
ance between the computed electronic exchange and correlation energies. In a simplified 
picture, the (negative) exchange energy is contributed by like spin electron pairs, while elec- 
tronic correlation arises from the interaction between electrons regardless of their spin. A 
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method which includes the exchange and neglects the correlation, as HF does, will favor 
high multiplicities by maximizing the number of electrons with the same spin. To the con- 
trary, experience shows that the combination of the exchange and the correlation terms in 
pure DFT pushes the balance toward low spin configurations.^^ Attempts to improve the 
spin state energetics description of density-functionals have mostly been based on hybrid 
Hartree-Fock/DFT schemes which combine the exchange of HF with the exchange 
and correlation obtained from DFT in proportions obeying empirical considerations. This 
approach, however, has given no universal functional capable to provide accurate splittings 
in every case. In general, those functionals offering a good description of the high spin 
species fail when tried out on low spin complexes, and vice versa.- Among them, B3LYP 
is seemingly the one with the best average performance up to now, yet exhibiting serious 
inaccuracies in the five coordinated models already discussed. 

In the present study, we propose the DFT+U approach as an alternative to the standard 
ab-initio techniques for a reliable description of the low lying states of iron heme complexes. 
The LDA+U or GGA+U method (more generally denoted as DFT+U) was originally de- 
signed within the density-functional theory framework for the treatment of strongly corre- 
lated materials.—"— Only very recently researchers have started to apply it to molecular, or 
mixed solid-molecular systems, with extremely promising results.—"— This approach cor- 
rects the tendency to overhybridize and delocalize electronic orbitals — ultimately originating 
in the presence of self-interactions in the exchange-correlation functionals — by introducing 
a term that penalizes fractional occupancies. We note in passing that in our present imple- 
mentation we explore the possibility of a [/ that is not a best-fit parameter, but an intrinsic, 
ab-initio linear-response property of the system chosen. However, this approach does not 
prove totally satisfactory, as in the case of low spin complexes it leads to values of U lying 
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1 or 2 eV above the optimal ones. We show that with the inclusion of a single parameter, 
DFT+U recovers the correct multiplicities of the five coordinated models where DFT and 
hybrid methodologies are in disagreement with more elaborated techniques or experimental 
data. Moreover, there is no impairment with respect to GGA functionals in those cases for 
which DFT displayed the right behavior. Calculations of ligand exchange thermodynamics, 
spin transitions, and other properties, point to GGA+U as an appealing tool to overcome 
the limitations entailed by the use of DFT in the description of bioinorganic complexes, at a 
computational expense much lower than demanded by highly-correlated quantum chemistry 
methods. 

II. METHODOLOGY 

A. General framework 

All calculations reported in this work have been performed with the public domain 
PWSCF and CP codes included in the Quantum-Espresso distribution,— based on density- 
functional theory, periodic-boundary conditions, plane-wave basis sets, and pseudopotentials 
to represent the ion-electron interactions. The PBE exchange-correlation functional^ has 
been used in combination with Vanderbilt ultrasoft pseudopotentials,— with the Kohn-Sham 
orbitals and charge density expanded in plane waves up to a kinetic energy cutoff of 25 and 
200 Ry respectively. 

B. The DFT-I-U approach 

The present implementation of DFT+U stems from the early contributions by Anisimov 
and others,—"— who proposed to correct the failures of the LDA functional in dealing with 
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the strongly localized d or f electrons of transition metal ions. An on-site correction was 
thus constructed to account for strong electronic correlations poorly described within the 
local-density or generalized-gradient approximations, and formulated as follows: 

Edft+uW)] = EDFTW)]+Eu[W^m']] = EdftW)]+ Ehu B[W^m']]- EDc[{n''']] (1) 

where p(r) is the electronic density, n^^/ are generalized atomic orbital occupations with 
spin a associated to the / atom, and rJ" is the sum of the occupations corresponding 
to all eigenstates, Y^mT^mm'- E£,FT[p{y)] is the standard LDA or GGA energy functional, 
and EHUB[{nl^,^i}] represents the "correct" on-site correlation energy. Since £'DFT[p(r)] 
already contains an approximate correlation contribution, a term intended to model such a 
contribution, i?£)c[{n^'^}], must be subtracted to avoid double counting. 

In this work, we resort to the rotationally invariant formulation of DFT+U introduced 
by Liechtenstein et alM- and later simplified by Dudarev and his coworkers ,"22. in which the 
non sphericity of the electronic interactions and the differences among the interactions in 
like-spin and unlike-spin channels are neglected. With these assumptions, the correction to 
the energy functional can be written 

EuiW^m'}] = ^ E E Km - E n'mm'nl:>J = ^ E Tr[^"^(l - ^'n] (2) 

where U is the Hubbard parameter describing on-site correlations. In principle, different 
definitions for the occupation matrix are possible, which in turn will determine different 
values for U. In this case we define 

ni:m' = j:MrM{<PL'm (3) 

u 

with /j, the weight of the electronic state z/, (pl^ the valence atomic orbital \lm) of atom 
/, and i/j^ the one electron wavefunction corresponding to the state u with spin a. The 
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diagonalization of the occupation matrices leads to the following expression for the energy 
correction: 

^^[{<^'}] = |EEAf'^(l-A^)- (4) 

1,(7 i 

Equation (4) clearly reflects the nature of the correction, which imposes a penalty (mediated 
by U) for fractional occupations, thus favoring either fully occupied or empty orbitals (A ~1 
and A ~0, respectively). We note that under this definition, U corresponds to the difference 
U — J a.s utilized by Anisimov and other researchers.—"— For example, the adoption of 
U=4 eV in the present calculations is comparable to a ?7 of 5 eV in combination with a 
J of 1 eV in the work of RoUmann.— Whereas in recent applications U is considered a 
fitting parameter,—"— here we obtain it from the spurious curvature of the DFT energy as 
a function of the occupations. As shown by Cococcioni and de Gironcoli,-- the value of U 
can be estimated as the difference between the screened and bare second derivative of the 
energy with respect to the occupations: 



U = r,,j., HTTTT^- (5) 



In particular, we are interested in the self-consistent [/, which we will call t/^c, originating 
from the curvature of the DFT+U ground state.— To compute Use-, a few linear response 
calculations must be performed at a finite Um, each one yielding a corresponding Uout- It 
can be shown that there is a linear dependence between [/j„ and Uout) from which Use can 
be extrapolated:— 

jT d Equad J J Uin 

'Jout — -TT, — TT^- — Use • yO) 

Equad groups all electronic terms within the DFT+U functional that have quadratic de- 
pendence on the occupations, whereas m can be interpreted as an effective degeneracy of 
the orbitals whose population is perturbed. This procedure, which allowed us to attain 



an improved description of the multiplet splittings and bonding in Fe dimers and FeO re- 
lated species,— is the one adopted here to calculate a self-consistent U parameter for the 
iron-porphyrin system. Another criterion has also been explored, requesting that a linear 
response calculation at a finite U returns this same value of U at the output, e.g. Uin = Uout- 
The parameter fulfilling this criterion will be hereafter denoted U'^^. This second criterion is 
not as appealing as the first one, since Use seems the "right definition" for self-consistency. 

III. RESULTS AND DISCUSSION 

In this section DFT-I-U results are presented on four heme complexes: Fe^^P(Im), 
Fe"^P(Cl), Fe"P(Im)(02), and Fe"P(CO). In the first two cases, DFT calculations fail to 
predict the high-spin nature of the system. The other two are examples of low-spin hemes 
whose electronic and geometrical properties are, in principle, correctly captured by standard 
density-functional simulations. These four case studies were chosen for their respective rel- 
evance in bioinorganic chemistry, and to illustrate the performance of the DFT+U method 
on heme models exhibiting a variety of coordination modes and multiplicities. 

A. Five coordinated heme-imidazole complex 

The Fe^^P(Im) system, depicted in Fig. [2], has been the target of numerous computational 
studies, inasmuch as it appears as the natural model to represent the unbound active site of 
several heme proteins, in particular hemoglobin and myoglobin As mentioned above, 
the ground state of this compound has been experimentally characterized as a quintuplet 
(S=2), whereas DFT calculations yield a triplet ground state (S=l). In Table II the energetic 
separations between the low lying spin states resulting from DFT and DFT+U are compared. 
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FIG. 2: Structure of the five coordinated Fe^^P(Im) complex. 

According to PBE, the triplet is around 8 kcal/mol more stable than the quintuplet or the 
singlet; similar gaps are obtained with the BP86 exchange-correlation functional.— On the 
other hand, the transition energy between the triplet and the quintuplet is reduced to nearly 
2 kcal/mol if B3LYP is used.— As previously noted, the introduction of the HF exchange in 
the DFT functional stabilizes high multiplicity states, but in this case B3LYP is still unable 
to provide the right splittings. 

The formalism summarized in equation (6) gives for this system a Use of 3.9 eV and a Ug^ 
of 2.5 eV (Fig. [3]). If any of these values are adopted, DFT+U restores the experimental 
ordering of the spin states. The total energies of the low lying states as a function of U are 
depicted in Fig. HI The increase of the U parameter equalizes the triplet and the quintuplet 
energies, producing a spin crossover at [/ ~ 2 eV. At higher values of U, the quintuplet 
remains the ground state. 

Fig. O highlights the effect of U on the electronic symmetry of the d states in the heme 
porphyrin. Spin occupations were computed by projecting the electronic wavefunctions on 
the atomic orbitals of the iron, as prescribed by equation (3). The upper panel of Fig. [5] 
represents the occupations of the minority spin manifold in the quintuplet state. Iron(II) is a 
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FIG. 3: Linear response calculation of the U parameter on the quintuplet state of the Fe P(Im) 
complex. 
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FIG. 4: Total energy of the low lying spin states of the Fe^^P(Im) complex as a function of the U 
parameter. 

ion and the quintuplet bears four unpaired electrons, therefore the sum of the occupations 
on the minority spin channel should be around 1. It is not exactly 1 because the eigenstates of 
the complex do not correspond to pure d atomic orbitals, but are instead strongly hybridized. 
Yet, it is possible to assign the electronic configuration of the system in terms of d atomic 
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FIG. 5: Occupations of the minority spin manifold in the quintuplet state of the Fe P(Im) complex. 
The lower panel shows the total energy of the two lowest accessible spin states as a function of U. 

orbitals, depending on whether the occupations are close to or 1. Note that by convention, 
the nitrogen atoms of the porphyrin ring are placed on the xy plane, with the x and y axes 
oriented along the Fe-N bonds. Pure DFT {U=0), and DFT+U with U < 2 eV, converge 
to the {dz^y {dxyY (dyr)^ {dx2^y2y state.— This is the same configuration as reported by 
Spiro and coworkers from B3LYP simulations.— The increase of U above 2 eV stabilizes the 
{dz^y {dxyY (dj,)^ {dx^-yiY state, which is the one experimentally assigned to Fe"P(Im).- 
Interestingly enough, this change in configuration is associated with an inversion in the 
relative energy of the triplet and the quintuplet, which now becomes the ground state. 

The examination of the optimized geometry at a finite U of 3.9 eV shows an agreement 
with the experimental data at least as good as pure DFT does. The structural parameters 
most affected by the U correction are those in the vicinity of the metal center, presented 
in Table III. The out-of-plane displacement of the iron, dp^^p, is the distance of the iron to 
the average plane defined by the four nitrogen atoms of the porphyrin ring. The interplay 
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FIG. 6: Distance of the iron to the average plane defined by the four nitrogen atoms of the porphyrin 
ring as a function of U in Fe^^P(Im). The shaded area encompasses the experimental region. 

between spin state and dpe-p, often involved in the dynamics of the heme protein — as for ex- 
ample in the allosteric mechanism of hemoglobin^i — has been characterized experimentally^ 
and theoretically.—*^ Table III contrast this and other optimized structural parameters with 
the experimental data available for the synthetic model compound Fe^^TPP(2-MeIm) (TPP: 
tetraphenyl porphine; 2-MeIm: 2-methyl imidazole).-- Fig. [6] shows the dependence of dpe-p 
on U, with the shaded part of the graph indicating the experimental region. 

In summary: the U term favors the stabilization of the {dxyY configuration — deemed the 
experimental ground state of the complex — rendering this (quintuplet) state the lowest in 
energy, as can be seen in the lower panel of Fig. [51 

B. Iron(III) porhyrinato chloride 

The low lying accessible electronic states for the penta-coordinate Fe^^^P(Cl) complex, 
with five d-electrons, are the sextet, the quartet, and the doublet (S=5/2, S=3/2 and S=l/2, 
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respectively). Similarly to the situation discussed in the previous section, DFT can not re- 
produce the high spin character of the system, which has been established experimentally.— 
Using a battery of ab-initio methods, Ghosh and collaborators have explored this complex 
in depth.-i^i^^ They found that, while the PW91 exchange-correlation functional yields a 
quartet state 8.1 kcal/mol more stable than the sextet, B3LYP provides nearly identical en- 
ergies for both configurations. Higher- level CASPT2 calculations and CCSD(T) simulations 
on a smaller model system are consistent with experiments, placing the sextet almost 20 
kcal/mol below the quartet.— These results are summarized in Table IV. 

Fig. [7] shows that, as seen in Fe^^P(Im), the U term stabilizes the highest multiplet in 
Fe^^^P(Cl). A spin inversion is verified at [/ ^ 1.5 eV, rendering the sextet as the ground 
state. A value of Use equal to 4.0 eV is obtained, which leads to a sextet-quartet transition 
energy of 9.2 kcal/mol. Table IV makes evident the poor performance of density-functionals 
to describe multiplet splittings in transition metals, capable of errors in the order of tens of 
kcal/mol. Despite its quantitative disagreement with the highly correlated methods (whose 
ultimate accuracy is, on the other hand, difficult to assess in this case), DFT-I-U succeeds 
in recovering the ordering of the spin states. 

C. Six coordinated oxyheme model 

The Fe^^P(Im)(02) system has been long identified as low spin in native proteins and 
in synthetic compounds.- Its importance as the oxygenated model of hemoglobin and myo- 
globin is reflected in the literature, which — aside from the experimental work — offers many 
computational studies addressing the electronic and structural aspects of the complex.— i^Sii^ 
The low spin nature of six coordinated iron porphyrins is in general correctly described by 
DFT, consequently with its trend to unstabilize high multiplicity states. In the particular 
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FIG. 7: Total energy of the lowest accessible spin states of the Fe^^P(Cl) complex as a function of 
the U parameter. 

case of Fe^^P(Im)(02), calculations with different functional, including B3LYP, indicate 
a singlet ground state of open-shell character. ^'^^ While the total spin of the molecule is 
zero, DFT calculations reveal partial spin densities localized on the d orbitals of Fe and the 
TT* orbitals of O2, integrating approximately to +1 and -1, corresponding to two unpaired 
electrons of opposite spin.-ii^ This open-shell singlet (o.s.s.) state can be interpreted as 
the result of an antiferromagnetic coupling between Fe^^P(Im) (S=2) and O2 (S=l), each 
retaining part of its magnetic character upon binding. 

DFT+U supports this picture: Fig. [S] depicts the spin density, Pspinij) = Pa(i") — P/3(r), 
computed at a finite ?7 of 4 eV. This figure is qualitatively equivalent to the one reported 
by Rovira and coworkers using pure DFT.— The impact of the U parameter on the total 
energies of the lowest accessible spin states is plotted in Fig. [9l The progressive increment 
of U further stabilizes the o.s.s. with respect to the closed-shell singlet and the triplet. On 
the other hand, the gap between the o.s.s. and the quintuplet is reduced, but the raise 
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FIG. 8: Spin density in Fe^^P(Im)(02) corresponding to an open-shell singlet, calculated with 
DFT+U. Lobes localized on the iron and on the O2 represent unpaired electron density of opposite 
spin. 

in U beyond 4 eV produces the dissociation of the Fe-0 bond in the o.s.s. before a spin 
crossing between these two states is observed. The effect of the on-site correction can also be 
examined through the absolute magnetization of the molecule, defined as / \pa{'r) — Pi3{r)\dr, 
a measure of the unpaired electron density in the system. Fig. illustrates how the on-site 
correlation affects the distribution of the unpaired electron density, reinforcing, in particular, 
the antiferromagnetic character of the o.s.s. The net effect of U is seemingly to thwart the 
coupling of the unpaired electrons of the molecular oxygen and the porphyrin, stabilizing 
the separate species, which is reflected in the elongation of the Fe-0 distance discussed 
below. The absolute magnetization augments from 1.8 at U=0 to 3.2 at U=A eV, a value 
placed halfway from that corresponding to the unbound system, comprised of a triplet plus 
a quintuplet (S=l+2). 

Extrapolation to the ?/-axis in the Uin — Uout plot (Fig. [H]) yields a Use of 5.9 eV, sensibly 
higher than in the previous examples. This value induces the rupture of the Fe-0 bond, and 
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FIG. 9: Total energy of the low lying spin states of Fe^^P(Im)(02) as a function of the U parameter. 
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FIG. 10: Absolute magnetization of the low lying spin states of Fe^^P(Im)(02) as a function of the 
U parameter. 

in consequence it is not possible to obtain a relaxed, bound complex associated with this Use- 
As seen in Fig. [T21 where the Fe-0 distance is plotted as a function of [/, the increase in the 
on-site correction provokes the elongation of the bond, eventually leading to the dissociation 
of the complex. The shaded area in the graph encloses the range of experimental Fe-0 
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lengths observed in different model compounds. It is important to emphasize that such 
values correspond to synthetic or natural compounds in which the O2 molecule is stabilized 
in the axial position by virtue of a second interaction on the distal side, namely an hydrogen 
bond or some kind of trapping or cage effect. In the absence of a distal cavity, oxygenation of 
iron (II) porphyrins under ordinary conditions has been rarely observed.— 1^ Pure DFT and 
B3LYP systematically overestimate the binding of O2, giving for the free heme energies in the 
range of 15-25 kcal/mol,— in direct contradiction with the experimental difficulty to isolate 
the oxygenated species. While part of this error is associated with the underestimation of 
the total energy of the quintuplet, such overbinding represents a major problem in the 
application of DFT to the calculation of affinity constants. Fig. [13] shows that DFT+U 
provides a more realistic oxygen affinity, the binding energy decreasing from 28 kcal/mol at 
U=0 to around 1 kcal/mol at ?7=4 eV. In the present case, Ui.^ turns out to be 1 eV lower 
than Use (Fig- [H])- A f/ of 5.8 eV, as obtained from equation (6), is too high to represent 
the thermodynamic and geometrical properties of the oxygenated complex. This is evincing 
a positive bias in the linear response approach, which will be manifested also in other low 
spin systems. We will come back to this issue later in the next sections. 

D. Five coordinated carboxyheme model 

As the last case study, we will briefly discuss the Fe^^P(CO) complex. Five and six 
coordinated carboxyhemes are low spin systems whose electronic and geometrical features 
are well reproduced by DFT, notwithstanding the overestimation of the CO binding energy, 
similarly to what is found with O2. The motivation to include the carboxylated complex in 
this study is therefore to assess the behavior of DFT+U in comparison with standard density- 
functional theory, in particular to examine if the former produces any detrimental bias in a 
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FIG. 11: Linear response calculation of the U parameter on the open-shell singlet state of the 
Fe^^P(Im)(02) complex. 
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FIG. 12: Fe-0 bond length as a function of the U parameter in Fe^^P(Im)(02). The shaded area 
encompasses the experimental region. 

case where the latter shows already a good performance. Table V contains computed and 
experimental values for a few selected properties of the carboxyheme. The general agreement 
between the simulations and the X-ray data is benefited from the U term, which not only 
provides a marginal improvement on the geometrical parameters of the complex, but also 
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FIG. 13: Energy of O2 binding to Fe"P(Im) as a function of the U parameter. 

corrects for the overbinding trend exhibited by pure DPT. At the same time, however, the 
enhancement of the on-site correlations closes the gap between the singlet and the quintuplet, 
to the extent that for U > 4 eV the latter becomes the most stable state (Fig. [13]). 

The response of the system to the change in orbital occupations resembles the case of 
Fe^^P(Im)(02), yielding for Use and U'^^ values of 7.2 and 5.3 eV respectively (Fig. [15]). Since 
neither experimental, nor reliable theoretical estimations of the spin transition energies are 
available, we can not evaluate the magnitude of the error in the self-consistent U obtained 
with each criteria. In principle, only values below 4 eV are consistent with the experimental 
singlet state, and so it is evident that both Use and U's^ suffer from some overestimation. 
Interestingly enough, the self-consistent U parameter calculated in a high spin configuration 
of the carboxyheme turns out to be significantly smaller, as depicted also in Fig. [13 On 
the other hand, linear response calculations on the low spin state of the Fe^^P(Im) system 
return values of Use above 6 eV (data not shown). This is indicating that the response of 
the system depends more on the multiplicity than on the particular geometry, coordination 
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FIG. 14: Total energy of the low lying spin states of Fe"P(CO) as a function of the U parameter. 
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FIG. 15: Linear response calculation of the U parameter on the singlet and quintuplet states of 
the Fe"P(CO) complex. 



mode, or the nature of the ligands. 
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IV. CONCLUSIONS AND FINAL REMARKS 



The inclusion of on-site correlations via a Hubbard term in DFT rectifies the trend of 
density-functionals to overstabilize low spin states in iron porphyrins. At variance with 
DFT, with Hartree-Fock, and with hybrid methods, which successfully describe some of 
the possible coordination modes of the complex but fail in the rest,- DFT+U is capable 
to provide the qualitatively correct splittings in low and high spin iron porphyrins at the 
same time, if the proper parameter is adopted. This improvement is also reflected in the 
geometry optimizations, and, more importantly, in more realistic binding energies to di- 
atomic ligands. The question of whether a hybrid functional with the proper exchange and 
correlation contributions would be capable to recover the spin state energetics of the full 
iron porphyrins series is still open. We have addressed this question in a previous article,- 
with no positive results. In our experience, the combination of the HF exchange with the 
GGA exchange-correlation leads to hybrid methods reflecting either the behavior of pure 
DFT — over stabilizing low spin states — or the behavior of Hartree-Fock — favoring high spin 
states. To the best of our knowledge, no hybrid has been reported that retains the best 
of both approaches in the description of iron porphyrins, but a more extensive search is 
probably needed before giving a definite answer. 

The application of the linear response calculation to low spin states leads to self-consistent 
Hubbard parameters 1 or 2 eV above the optimal ones. The linear response of the system 
appears to be more dependent on the spin state than on the coordination number or the 
identity of the hgands. In fact, plots of Uin versus Uout belonging to different complexes 
in the same spin state exhibit a high similarity. The reason for the overestimation of the 
Hubbard U in low spin configurations is not evident. Different extensions to the linear 
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response approach were explored, including the partition of U into U°' and to discriminate 
between both spin channels, and even between the five d states. Additionally, a J term to 
represent separately the on-site exchange was implemented. The modifications described 
above, however, produced little or no effect on the resulting Ugc- The investigation of other 
factors which could be responsible for these biases is in progress. In any case, values of U 
of 4 eV or slightly lower seem the optimal to reproduce the electronic, thermodynamic and 
structural properties of the heme compounds. 
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TABLE I: Experimental and calculated electronic ground states of five and six coordinated iron 
porphines (FeP), with the following axial ligands: O2, CO, imidazole (Im), and chloride. 



Six coordinated Five coordinated 





FeP(Im)(02) 


FeP(Im)(CO) 


FeP(CO) 


FeP(Im) 


FeP (CI) 


Experimental 


singlet 


singlet 


singlet 


quintuplet 


sextet 


Hartree-Fock 


quintuplet 


quintuplet 


quintuplet 


quintuplet 


sextet 


DFT-GGA 


singlet 


singlet 


singlet 


triplet 


quartet 


B3LYP 


singlet 


singlet 


singlet 


triplet 


quartet / sextet 



TABLE IL Spin transition energies (kcal/mol) for the low lying spin states of Fe P(Im) calculated 
with several density-functionals and with DFT+U, using [/sc=3.9 eV. 



Singlet Triplet Quintuplet 

PBE*^ 7.8 0.0 7.9 

DFT BP86'' 8.3 0.0 6.5 

BSLYP'^ 5.8 0.0 1.9 

DFT+U 20.9 4.9 0.0 



"Pseudopotential calculations with plane wave basis set. ^Pseudopotential calculations with plane 
wave basis set, from ref.^^. ^Gaussian calculations (VTZ basis), from ref.— . 
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TABLE III: Selected experimental and optimized structural parameters (A) in Fe^^P(Im). 



dpe-p Fe-Nporp/i Fe-Njmi 



Experimental'* 0.42 2.09 2.16 

DFT (PBE) 0.29 2.07 2.13 

DFT+U {Usc=3.9 eV) 0.43 2.11 2.19 



Data for Fe"TPP(2Me-Im) from ref.^. 



TABLE IV: Spin transition energies (kcal/mol) for the low lying spin states of Fe P(C1) calculated 
with highly correlated methods and density functional theory, including DFT+U (C/sc=4.0 eV). 



Quartet Sextet 
CASPT2" 19.6 0.0 
RCCSD(T)'' 16.1 0.0 
DFT-PBE 0.0 5.6 
DFT-PW91'^ 0.0 8.1 
DFT+U 9.2 0.0 



Ref.— . ''Calculations on a simplified model, ref/'^-. "^Ref.— . 
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TABLE V: Selected bond distances (A) and CO binding energy (kcal/mol) for the Fe"P(CO) 
complex calculated with DFT+U. 



U (eV) Fe-C(CO) Fe-Nporph C-0 AE 

0.0 1.69 1.99 1.17 -45.3 

3.0 1.71 2.00 1.16 -26.1 

5.0 1.74 2.01 1.16 -12.0 

Experimental" 1.77 2.02 1.12 



"The experimental model is axially coordinated to pyridine, ref.=. 
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